library(sf)
library(tidyverse)
library(tmap)
library(raster)
library(mapedit)
library(magrittr)
library(stars)
library(units)Unidad 3: Álgebra de Mapas y Geoprocesamiento
Manipulando datos espaciales
- La tabla de atributos
- Herramientas de cálculo
- Funciones y variables
- Trabajo con las geometrías
- Conversión de geometrías
- Conversión de modelos:
raster <-> vector
Paquetes utilizados en esta lección
Actividad 1. Clases, atributos y valores
- Cargue las capas ocurrencias y Mask del fichero data/gpkg/tapir.gpkg
tapir_mask <- st_read("data/gpkg/tapir.gpkg", layer = "mask")
tapir_ocu <- st_read("data/gpkg/tapir.gpkg", layer = "ocurrencias")- Explore los datos y determine qué tipo de dato es cada campo.
glimpse(tapir_mask)Rows: 1
Columns: 2
$ Name <chr> NA
$ geom <MULTIPOLYGON [°]> MULTIPOLYGON (((-81.47795 -...
Para el caso de
mask, solo tiene un campo y es tipocharacter
glimpse(tapir_ocu)Rows: 155
Columns: 7
$ Especie <chr> "Tpinchaque", "Tpinchaque", "Tpinchaque", "Tpinchaque", "Tp…
$ Confirmado <lgl> TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE,…
$ timestamp <dttm> 2008-08-18 09:48:42, 2009-09-29 19:37:12, 2009-02-12 19:47…
$ Longitud <chr> "78º 21' 24.12\" W", "77º 34' 24.13\" W", "77º 32' 54.13\" …
$ Latitud <chr> "0º 40' 38.24\" N", "2º 15' 8.24\" N", "2º 11' 38.26\" N", …
$ zona_UTM <chr> "18N", "18N", "18N", "18N", "18N", "18N", "18N", "18N", "18…
$ geom <POINT [°]> POINT (-77.6433 0.67729), POINT (-76.42663 2.25229), …
La capa de
ocurrenciaspresenta 6 campos y en su mayoría soncharacters, también haylogicalydatetime.
- ¿Que sistema de coordenadas tienen las capas?
La capa de mask, tiene un epsg: 4326, mientras que la capa ocurrencias, tiene un epsg: 4326.
- ¿Qué tipos de operaciones se pueden hacer con esos tipos de datos?
Con la columna tipo
logicalse puede realizar operaciones booleanas, somo sumar los verdaderos. Con las fechas probablemente una diferencia de fechas o conteo por un período de tiempo (e.g. por meses), una operación similar se podría hacer por Especie.
Actividad 2. Cálculos avanzados: Funciones personalizadas
Usando las capas del fichero “data/gpkg/muestreo.gpkg”:
muestreo <- st_read("data/gpkg/muestreo.gpkg")- Agregar una columna que contenga las coordenadas X y Y en el SRC UTM-17S.
muestreo$X <- st_coordinates(muestreo)[, "X"]
muestreo$Y <- st_coordinates(muestreo)[, "Y"]
knitr::kable(head(muestreo))| id | ta_media | Altitud | Forest_P_2010 | Forest_P_2000 | geom | X | Y |
|---|---|---|---|---|---|---|---|
| 0 | 11.6 | 3513.165 | 22.5325 | 18.28994 | POINT (763860.6 9677462) | 763860.6 | 9677462 |
| 1 | 18.8 | 1899.906 | 43.8950 | 25.13772 | POINT (746943.6 9598100) | 746943.6 | 9598100 |
| 2 | 18.8 | 1937.868 | 26.9850 | 9.45509 | POINT (749704.4 9644103) | 749704.4 | 9644103 |
| 3 | 16.9 | 2144.949 | 37.0575 | 24.38037 | POINT (732393.7 9612334) | 732393.7 | 9612334 |
| 4 | 11.8 | 3073.430 | 37.7225 | 38.11243 | POINT (675170.9 9598945) | 675170.9 | 9598945 |
| 5 | 15.2 | 2740.981 | 36.6025 | 20.14110 | POINT (779156.3 9688753) | 779156.3 | 9688753 |
- Agregar una columna nueva que contenga los pesos ficticios entre 70 y 120 kg, asignados de manera aleatoria a cada registro
muestreo$Weight <- sample(70:120, size = nrow(muestreo), replace = T)
knitr::kable(head(muestreo))| id | ta_media | Altitud | Forest_P_2010 | Forest_P_2000 | geom | X | Y | Weight |
|---|---|---|---|---|---|---|---|---|
| 0 | 11.6 | 3513.165 | 22.5325 | 18.28994 | POINT (763860.6 9677462) | 763860.6 | 9677462 | 74 |
| 1 | 18.8 | 1899.906 | 43.8950 | 25.13772 | POINT (746943.6 9598100) | 746943.6 | 9598100 | 96 |
| 2 | 18.8 | 1937.868 | 26.9850 | 9.45509 | POINT (749704.4 9644103) | 749704.4 | 9644103 | 109 |
| 3 | 16.9 | 2144.949 | 37.0575 | 24.38037 | POINT (732393.7 9612334) | 732393.7 | 9612334 | 99 |
| 4 | 11.8 | 3073.430 | 37.7225 | 38.11243 | POINT (675170.9 9598945) | 675170.9 | 9598945 | 104 |
| 5 | 15.2 | 2740.981 | 36.6025 | 20.14110 | POINT (779156.3 9688753) | 779156.3 | 9688753 | 113 |
- Agregar una columna nueva que contenga la zona UTM a la que pertenecen
GetUtmZone = function(x){
centroid = st_centroid(x)
longitude = st_coordinates(centroid)[,1]
latitude = st_coordinates(centroid)[,2]
zone_number = floor(((longitude + 180) / 6) %% 60) + 1
zone_letter = ifelse(latitude >= 0,
'N', #if true
'S') #else
return(paste0(zone_number, zone_letter))
}
muestreo$Zone <- st_transform(muestreo, 4326) |> st_geometry() |> GetUtmZone()
knitr::kable(head(muestreo))| id | ta_media | Altitud | Forest_P_2010 | Forest_P_2000 | geom | X | Y | Weight | Zone |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 11.6 | 3513.165 | 22.5325 | 18.28994 | POINT (763860.6 9677462) | 763860.6 | 9677462 | 74 | 17S |
| 1 | 18.8 | 1899.906 | 43.8950 | 25.13772 | POINT (746943.6 9598100) | 746943.6 | 9598100 | 96 | 17S |
| 2 | 18.8 | 1937.868 | 26.9850 | 9.45509 | POINT (749704.4 9644103) | 749704.4 | 9644103 | 109 | 17S |
| 3 | 16.9 | 2144.949 | 37.0575 | 24.38037 | POINT (732393.7 9612334) | 732393.7 | 9612334 | 99 | 17S |
| 4 | 11.8 | 3073.430 | 37.7225 | 38.11243 | POINT (675170.9 9598945) | 675170.9 | 9598945 | 104 | 17S |
| 5 | 15.2 | 2740.981 | 36.6025 | 20.14110 | POINT (779156.3 9688753) | 779156.3 | 9688753 | 113 | 17S |
- Agregue una columna nueva que contenga un identificador para los puntos que pertenecen al hemisferio sur y que caigan dentro del área que cubre la Demarcación hidrográfica Santiago
DHS <- st_read(
"data/gpkg/muestreo.gpkg", layer = "Demarcacion hidrogafica Santiago",
quiet = T
)
muestreo$dhs_sur <- st_intersects(muestreo, DHS, sparse = F) & muestreo$Zone == "17S"
knitr::kable(head(muestreo))| id | ta_media | Altitud | Forest_P_2010 | Forest_P_2000 | geom | X | Y | Weight | Zone | dhs_sur |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 11.6 | 3513.165 | 22.5325 | 18.28994 | POINT (763860.6 9677462) | 763860.6 | 9677462 | 74 | 17S | TRUE |
| 1 | 18.8 | 1899.906 | 43.8950 | 25.13772 | POINT (746943.6 9598100) | 746943.6 | 9598100 | 96 | 17S | TRUE |
| 2 | 18.8 | 1937.868 | 26.9850 | 9.45509 | POINT (749704.4 9644103) | 749704.4 | 9644103 | 109 | 17S | TRUE |
| 3 | 16.9 | 2144.949 | 37.0575 | 24.38037 | POINT (732393.7 9612334) | 732393.7 | 9612334 | 99 | 17S | TRUE |
| 4 | 11.8 | 3073.430 | 37.7225 | 38.11243 | POINT (675170.9 9598945) | 675170.9 | 9598945 | 104 | 17S | FALSE |
| 5 | 15.2 | 2740.981 | 36.6025 | 20.14110 | POINT (779156.3 9688753) | 779156.3 | 9688753 | 113 | 17S | TRUE |
Actividad 3. Edición de capas de líneas y polígonos
- Cargue la librería
mapedit
library(mapedit)Ejecute la función
mapedit::drawFeatures()Acérquese a la ciudad de Cuenca, específicamente al Parque Calderón (Parque central de Cuenca)
Trace puntos de los 2 o 3 edificios representativos que que pueda reconocerlos fácilmente alrededor del Parque.
Trace un rectángulo que cubra hasta tres cuadras alrededor del parque calderón.
Trace las vías hasta dos cuadras alrededor del Parque.
Guarde el resultado en una capa, con SRC EPSG:4326 (Latitud y Longitud).
parque_calderon <- drawFeatures()
knitr::kable(head(parque_calderon))| X_leaflet_id | feature_type | geom |
|---|---|---|
| 1622 | marker | POINT (-79.00394 -2.898353) |
| 1628 | marker | POINT (-79.0053 -2.897426) |
| 1634 | marker | POINT (-79.00459 -2.898053) |
| 2017 | rectangle | POLYGON ((-79.00741 -2.9017… |
| 2058 | polyline | LINESTRING (-79.0069 -2.897… |
| 2082 | polyline | LINESTRING (-79.00678 -2.89… |
Agregue una columna con el nombre o referencia de las vías dibujadas.
Explore el resultado y discuta sobre las propiedades de esta capa.
Para esto visualizaremos el mapa con tmap y señalaremos las ids de las vías para asignar el nombre que le corresponde:
tmap_mode("view")
filter(parque_calderon, st_geometry_type(parque_calderon) == "LINESTRING") |>
qtm()Asignación de nombres a las vías
parque_calderon %<>%
mutate(
vias_id = case_when(
`X_leaflet_id` == 2358 ~ "Gran Colombia",
`X_leaflet_id` == 2082 ~ "Simón Bolívar",
`X_leaflet_id` == 2058 ~ "Mariscal Sucre",
`X_leaflet_id` == 2154 ~ "Presidente Córdova",
`X_leaflet_id` == 2174 ~ "Juan Jaramillo",
`X_leaflet_id` == 2238 ~ "Padre Aguirre",
`X_leaflet_id` == 2266 ~ "Benigno Malo",
`X_leaflet_id` == 2292 ~ "Luis Cordero",
`X_leaflet_id` == 2316 ~ "Presidente Borrero",
T ~ "Otras entidades"
)
)Es un tipo de multi-geometría, es decir tiene una mezcla de polígonos, puntos y líneas. A pesar de que R lo puede procesar, otros programas como Qgis pueden necesitar que estás geometrías estén separadas.
- Separe las entidades por geometrías y guarde como capas separadas en un fichero
*.gpkg.
filter(parque_calderon, st_geometry_type(geometry) == "LINESTRING") |>
st_write("data/gpkg/parque_calderon.gpkg", layer = "pc_lines")
filter(parque_calderon, st_geometry_type(geometry) == "POINT") |>
st_write("data/gpkg/parque_calderon.gpkg", layer = "pc_points")
filter(parque_calderon, st_geometry_type(geometry) == "POLYGON") |>
st_write("data/gpkg/parque_calderon.gpkg", layer = "pc_polygons")Actividad 4. Geometrías derivadas
- Cargue la capa puntos disponible en el fichero quinta_balzay.gpkg
quinta_balzay <- st_read("data/gpkg/quinta_balzay.gpkg")- Explore las propiedades de la capa y haga las preguntas necesarias para entender los atributos.
- Genere una o varias capas de según corresponda. Si es necesario convierta esas geometrías a geometrías más complejas.
Para evitar que las geometrías se agrupen incorrectamente, hay que aplicar el argumento
do_union = Fensummarise()
Formación de geometrías a polígonos
polygons_qb <- quinta_balzay |>
filter(feature_type == "polygon") |>
group_by(id) |>
summarise(do_union = F) |>
st_cast("POLYGON")Formación de geometrías a líneas
lines_qb <- quinta_balzay |>
filter(feature_type == "polyline") |>
group_by(id) |>
summarise(do_union = F) |>
st_cast("LINESTRING")Extracción de puntos
points_qb <- quinta_balzay |>
filter(feature_type == "point")qtm(polygons_qb, fill = "id") +
qtm(lines_qb) + qtm(points_qb)- Explore los resultados y comente lo que descubrió.
Los polígonos hacen referencia a los aularios y a las oficinas cerca de hídrica, las líneas a las vías dentros del campus balzay y los puntos a zonas importantes de la casa antigua de la quinta balzay (e.g. laboratorios).
Actividad 5. Raster a vector
- Cargue la capa de Temperatura_3.tif y convierta las clases a polígonos
- A partir de la capa Temperatura_2.tif, genere isolíneas cada 5 grados de temperatura.
at_3 <- raster("data/tif/Temperatura_3.tif")
at_2 <- raster("data/tif/Temperatura_2.tif")A continuación se recorta las capa de temperatura con la de muestreo, para ahorrar tiempo en el procesamiento de las capas.
- Isolíneas cada 5 grados
at_2_iso <- at_2 |> crop(muestreo) |>
cut(breaks = seq(-5, 35, 5), include.lowest = TRUE, right = FALSE) |>
raster::rasterToContour() |>
st_as_sf()- Polígonos de las clases
at_3_poly <- at_3 |> crop(muestreo) |> st_as_stars() |> st_as_sf(merge = T)qtm(at_2_iso, lines.col = "level")
qtm(at_3_poly, fill = "Temperatura_3")
qtm(at_3_poly, fill = "Temperatura_3") +
qtm(at_2_iso, lines.col = "level")Geoprocesamiento y Álgebra de mapas
- Geoprocesamiento (Vector)
- Análisis de solape
- Análisis de proximidad
- Análisis de agregación
- Álgebra de mapas (Raster)
- Operaciones locales
- Operaciones zonales
- Operaciones globales
Actividad 5. Raster a vector
Una laguna de interés de un estudio se ecuentra en la siguiente coordenada: POINT(-79.27317 -2.85857). Usando las capas del fichero geoprocesos.gpkg, encuentre las respuestas a las siguientes preguntas:
Creamos el punto espacial e importamos la capa de geoprocesos.gpkg
laguna_point <- "POINT(-79.27317 -2.85857)" |> st_as_sfc(crs = 4326) |>
st_as_sf() |> st_transform(crs = 32717)
water_bodies <- st_read("data/gpkg/geoprocesos.gpkg", layer = "cuerpos de agua") |>
st_transform(crs = 32717)- ¿En un radio de 6 Km desde el punto, cuántos cuerpos de agua naturales existen?
st_is_within_distance(
x = laguna_point, y = water_bodies,
dist = units::set_units(6,"km"),
sparse=FALSE
) %>% sum()[1] 126
Existen 126 cuerpoos de agua al rededor de los 6km de la laguna de interés.
- ¿Cuál es el área que cubren todos ellas?
lag_buff <- st_buffer(laguna_point, dist = units::set_units(6, "km"))
lag_buff |>
st_area() |> units::set_units("km^2")113.0457 [km^2]
Tiene un área total de 113.05 km\(^2\)
- ¿Cuál es el rango de temperaturas que existen en la zona de influencia? (Use la capa de Puntos de muestreo disponible en muestreo.gpkg)
qtm(muestreo) +
qtm(lag_buff)muestreo[lag_buff, ] |> pull(ta_media) |> mean()[1] 8
Tiene una media de 8°C
- ¿Cuánto espacio del radio de 6 km, no está cubierto por cuerpos de agua?
st_difference(
lag_buff, st_union(water_bodies)
) %>% st_area() %>%
units::set_units("km^2")106.433 [km^2]
No están cubiertos por agua 106.4 km\(^2\)
Actividad 6. Álgebra de mapas
En un proyecto nuevo ya sea en QGIS o R cargue la capa temperatura_2.tif. Usando las herramientas necesarias genere una capa como la que se muestra en la imagen.

tm_shape(at_2) +
tm_raster(
palette = "-RdYlBu",
midpoint = 14,
style = "cont"
)at_levs <- cut(
at_2, breaks = c(-Inf, -1, 5, 15, 22, 25, 50)
)
tm_shape(at_levs) +
tm_raster(palette = "viridis", midpoint = 3)Actividad 7. Álgebra de mapas 2
**Usando los recortes de Landsat 8 (clip_RT_LC08_L1TP_010063*), calcule el NDVI para toda la zona de estudio. El NDVI se calcula mendiante la ecuación:**
\[ NDVI = \frac{NIR - Red}{NIR + Red} \]
Importamos los datos y necesitamos las bandas 5 y 4
lc_f <- list.files(
"data/tif", pattern = "clip_RT",
full.names = T
)Cálculo de NDVI
ndvi_fun <- function(nir, red) {
(nir - red) / (nir + red)
}
ndvi <- ndvi_fun(
nir = raster(lc_f[4]), red = raster(lc_f[3])
)